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Constraining corotation from shocks in tightly-wound 
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ABSTRACT 

We present a new method for estimating the corotation radius in tightly wound spiral 
galaxies, through analysis of the radial variation of the offset between arms traced by 
the potential (P-arms) and those traced by dust (D-arms). We have verified the predic- 
tions of semi-analytical theory through hydrodynamical simulations and have exam- 
ined the uniqueness of the galactic parameters that can be deduced by this method. 
We find that if the range of angular offsets measured at different radii in a galaxy 
is greater than around 7r/4, it is possible to locate the radius of corotation to within 
~ 25%. We argue that the relative location of the P- and D-arms provides more robust 
constraints on the galactic parameters than can be inferred from regions of enhanced 
star formation (SF-arms), since interpretation of the latter involves uncertainties due 
to reddening and the assumed star formation law. We thus stress the importance of 
K-band studies of spiral galaxies. 
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1 INTRODUCTION 

Theories of the structure, origin, development and effects of 
spiral arms in galaxies have occupied researchers ever since 
such structure was discovered in 1845. Many galaxies are 
observed with tightly-wound, regular spiral arms covering a 
significant fraction of their disc. The commonest model of 
this pheno menon is based on the L in-Shu hypothesis of spiral 
structure iLin. Yuan fc Shulll969l) . which states that quasi- 
steady spiral patterns exist in the stellar disc, rotate at a 
pattern speed flp, and persist for many orbits. The orbits of 
the stars are organised into 'kinematic density waves', which 
increase the stellar density in a spiral pattern. The stars 
themselves drift through the pattern inside the corotation 
radius, as they orbit with an angular velocity Q > f2p. 

The spiral density wave picture is by no means the only 
theory of spiral structure, but provides a very successful 
description. The possible origins, driving mechanisms and 
long-term support of such density waves has been the sub- 
ject of a great deal of literature, and thi s subject is beyond 
the sc ope of this paper; s ee reviews bv iKaplan fc Pikelnerl 
il974h and lToomrej ( 1977). It is generally accepted that spi- 
ral waves, if they do not naturally arise through instabilities, 
may be driven by nearby companions or central bars, and 
that under the right conditions they may persist for many 
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rotations, and it is from this starting point that we proceed 
in this paper. 

Galaxies that are optically classified as spirals usually 
have a clear spiral structure with one dominating mode, 
and arms are usually l ogarithmic in shape (e.g. they have a 
constant pitch angle 1 l| Kennicutt||l98l|: [Kennicutt fc H oded 
ll982tlGarci'a Gomez fc Athanassoulall 19931 lMall2002lh How- 
ever, galaxies not classified as spirals are now known com- 
monly to show spiral structure in their stellar discs. The old 
stellar population, which will trace the mass of the disc, can 
be observed in lo nger wavelengths such as the K or K' bands. 
ISchweizerl Jl976l) found 'broad spiral patterns' in the discs of 
six spiral galaxies, concluding tha t smooth spiral structure 
exists in their mass distributions. lElmegreen fc Elmegreenl 
il984h examined several galaxies of different types, sug- 
gesting a class of galaxy which has flocculent structure at 
blue wavelengths a nd a spiral structure in K. More recently, 
iBlock et alJ <ll994l) showed that optical classification does 
not constrain the structure of stellar discs, finding smooth 
spiral structure in galactic discs independently of their op- 
tical structure. Many galaxies optically classified as floc- 
culent, when observed in K' , show regular spiral structur e 
<ThornlevHl993 iGrosbfll fc Patsislll998l ISeiear et alJl2003h . 
Most of these have a central bar ] and s ome also have a neigh- 
bour iSeiear. Chornev fc Jameal2003h . The presence of reg- 
ular spiral structure, but flocculent structure in blu er wave- 
lengt hs, has also arisen in numerical simulations llBermanl 
2002). We make the distinction between three separate arm 
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components: the peaks of the underlying stellar mass distri- 
bution, approximately tracing the perturbations in the grav- 
itational Potential (the P-arms), the shock fronts of the gas, 
which will generally be observed from Dust (the D-arms), 
and the ridges of enhanced Star Formation (the SF-arms). 

Many tightly-wound spiral galaxies are known. Assum- 
ing that the model of a rotating spiral potential is accu- 
rate, the parameters of that potential are fundamental to 
the galaxies themselves. An understanding of the structure 
of spiral arms is important in considering the processes af- 
fecting the ISM flowing through them. In this paper, we 
focus on the determination of the corotation radius. 

We discuss the formation of single-shock flows in the 
gaseous disc, and show how the angular offset between these 
shocks (the D-arms) and the minima of the gravitational 
potential (the P-arms) could be used to constrain corota- 
tion. In section |3 we introduce the model spiral potential, 
to which the gas responds. Sectional investigates the result- 
ing gas flow using a well known semi-analytical approach, 
and compares the results to those obtained by numerical 
methods. In section 2] different components of spiral arms 
are discussed, and the offset function O is introduced. Vari- 
ous ways to constrain corotation are summarised in section 
and the use of the offset function to provide an additional 
constraint is explained. Finally, our findings are summarised 
in section HJ 



2 MODEL SPIRAL POTENTIAL 

The rigidly rotating potential of a spiral galaxy can be de- 
composed into a static, axisymmetric part and a spiral per- 
turbation. We define the gravitational potential V (in galac- 
tocentric polar coordinates R, 9) as 

V{R,9,t) =V R (R) + V s {R,e,t) 

where Vr is the axisymmetric potential and Vs is the spiral 
perturbation. 

The axisymmetric part represents the unperturbed ve- 
locity curve. In this paper, the following velocity curve is 
used: 

v(R) = Vmaxyt' F b e b Rexp(— e b R) + 1 — exp(— e d R) 

where t; ma x is the limit of the circular velocity at large 
radius, ed and e b are the inverse disc and bulge scale 
lengt hs (respectively) and F b is th e bulge strength param- 
eter dContopoulos fc Gro sbol 1986). The potential Vr that 
produces this velocity curve is 

V n (R) = «Lx (inR + E^R) - F h e~^ R ) 

where the exponential integral Ei is defined as 

, . dit 

E 1 (x) = 

A spiral perturbation with one mode is given by the 
general form 

Vs = A(R,t) cos(x) 

where the spiral phase \ is defined as 

X = HR) - m{6 - flpt) 
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0, 
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Figure 1. Potential surface with parameters F b — u , c d 
kpc, e^ 1 = 10 kpc, sini = 0.2, m = 2 and Fq = .6 at Rq = 
8.5 kpc. A very high value of F is chosen to produce a spiral 
perturbation visible on this figure. 



and the pattern has m arms. The pattern rotates rigidly 
at angular velocity Qp. The radial shape of the arms is set 
by the function &(R), and at constant radius, the potential 
is sinusoidal in 6 with period 2n/m. The overall amplitude 
is set by A(R,t). The equipotentials of each mode are spi- 
rals, which have an inclination i(R) to circles. Trailing loga- 
rithmic arms, with a constant inclination i, correspond the 
choice 
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It is convenient to define the perturbation strength 
F(R) as the ratio of the amplitude of the perturbation force 
to the axisymmetric force: 



F 
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where the circular velocity at radius R is v. The radial and 
time variation of the amplitude can be chosen arbitrarily. A 
suitable radial choice is to specify llContopoulos fc Grosbe)l| 

A(R,t) = A Rexp(-e B R) 

where e s is the inverse spiral scale length, and the over- 
all amplitude at any time is set by Aq. In this paper we 
are interested in steady solutions in which the potential has 
reached a stable state, so that A is independent of time. Aq 
can most conveniently be set by requiring that the maxi- 
mum relative spiral strength takes a specified value Fq at 
some radius Ro- 

An example potential surface, with an exaggerated spi- 
ral perturbation, is shown in figure^ 



3 RESPONSE OF GASEOUS DISC 

A relatively small perturbation can have a significant effect 
on the flow of gas through it. Large deviations from the 
circular flow can lead to the formation of shocks. In this 
section we examine the formation of single-shock solutions 



in the flow of isothermal gas through the spiral potential. 
Since the intent is to examine the large-scale behaviour, and 
ignore small-scale processes, approximations of initial uni- 
formity and isothermality are appropriate. Here, the sound 
speed used is an 'effective' sound speed equal to the velocity 
dispersion of the ISM, a g ood approxi mation for the time 
and length scales involved (ICowielll98d) . 

3.1 Semi-analytical technique 

In this paper we re peat the procedure used in 
IShu. Milione fc Roberts) il973li for finding non-linear 
solutions (containing a shock) to the asymptotic equations 
of gas flow under a spiral perturbation. The derivation of 
th e equations an d de scriptions of the procedure are given 
in lRobertsI Jl969ft and lShu et all (I1972T) and are only briefly 
described here. 

Solutions are expressed in 'curvilinear' coordinates 
(?7,£) that vary perpendicularly to and parallel to the 
equipotentials, respectively (see figure |5J. These are defined 
in a stationary frame of the spiral perturbation. They are 
defined by 1 

drj = —kdR + md9 
riff 

dt = -m-?-- kRd6 
R 

where 

_ d${R) _ -m 
dR Rtani 

where the last result is for the case of logarithmic arms. A 
circle at constant R will therefore pass through m periods 
of 77. Specifically we define 

T) = -X + * (!) 

so that 7} = always corresponds to the potential minimum. 
The passage from one arm to the next corresponds to a 
period of 2ir in 77, regardless of the number of arms. 

The velocity of the gas is written in components (u v , tij ) 
in this coordinate system. These are then expressed as the 
sum of the 'base' flow (w n o,W{o) (the equilibrium circular 
flow when the perturbation is not present, in the station- 
ary frame of the perturbation) and a perturbed velocity 
(u v i,U(-i). In the specific case of logarithmic arms the co- 
ordinates can be defined as follows: 
71% 

Tj = ; ln(e a R) + m(6 - Q P t) + tv 

tan 1 

£ = — mln(e s _R) + mtani(# — Qpt). 

The differential equations in u v and u^, and the tech- 
nique for locating appropriate solutions, are given in ap- 
pendix^ The boundary condition for any solution is that 
it be periodic in 77 with period 2n. Solutions are found along 
streamlines under the approximation that they are at nearly 
constant radius, and the surface density a is then found from 
the velocity. The solution at a given radius is controlled by 
seven parameters, under the approximation that they are 
constant along a streamline: 

T hese coordinate s differ slightly to those in IShu et al ] <1972fl 
and lRobertsl <1969l) . 
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Figure 2. The coordinate system used in this paper. 



m (the number of arms) 

Qp (the pattern speed of the arms) 

i (the local pitch angle) 

F (the local relative arm strength) 

v (the unperturbed circular velocity) 

k (the local epicyclic frequency) 

a (the effective sound speed of the gas) 

Base-supersonic flow is defined as any flow with u v o > a, 
so that in the stationary frame of the arms, the unperturbed 
flow perpendicular to the equipotentials is supersonic. Sim- 
ilary, base-subsonic flow has u v o < a- In general, if the flow 
in a galaxy is mostly base-supersonic, there will be a region 
(the base-subsonic region) in which it is not. This, of course, 
occurs either side of the corotation radius in a band where 
the relative velocity R(Q. — fip) is sufficiently small. In this 
region, the response of the gas, and the asymptotic flow so- 
lutions, will be qualitatively different. The base-subsonic re- 
gion is bounded by the radii satisfying i?sini(Sl — f2p) = ±0. 

If for the whole solution u v < a (entirely subsonic) or 
u v > a (entirely supersonic), the solution will not contain 
a shock. This will occur under conditions where the sound 
speed a is very high or low, or the strength F of the pertur- 
bation is sufficiently small. In the limit of small F, the solu- 
tions reach the first-order response, which is sinusoidal in u 
and v (and hence a), such that u n oc — cos 77 and oc sin 77. 
The maximum density (and minimum of u v ) occur at the 
centre of the potential well, i.e. at 77 = 0. 

If the strength F of the perturbation is gradually in- 
creased, starting from an entirely subsonic solution, then at 
some point the flow will be sufficiently perturbed that u v will 
pass the sound speed a. Any solution must now contain a 
sonic point, at which the gas accelerates to supersonic veloc- 
ities, and a shock. Such solutions are no longer symmetric. A 
consistent solution consists of a closed curve in the (u^,u v ) 
plane containing a single isothermal shock, with a periodic- 
ity in 77 of 2-k. Such a solution only exists for certain choices 
of parameters. 

Figure|3]shows an example velocity curve in (u^, u v ) for 
a solution in a tightly wound spiral, and the corresponding 
density profile through the shock is shown in figure [1] 
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Figure 3. Velocity curve in [uc, u,,) plane of the solution at Rq = 
8 kpc and F = 0.05. The cross marks the unperturbed velocity. 
The sonic point is marked by a circle. The shock is represented 
by the dashed line. 




V 



Figure 4. Surface density profile cr(r]) corresponding to the so- 
lution shown in figure |3 The sonic point is marked by a circle. 



Table 1. Parameter values of the standard model. 








-i 


1.5 kpc 




10 kpc 


vrt 


2 


sin i 


0.1 


tt P 


13 km s 


a 


8 km s 


F(Ro) = F 


0.05 


v(Ro) = vq 


220 km s 


Ro 


8.5 kpc 



only if F is not much larger than sini. In our experience, 
solutions are difficult or impossible to find for F>0.1 and 
sin i> 0.2 (with all other parameters as indicated in table 

m ~ 

Within parameters that satisfy these approximations, 
solutions still may not exist. In particular, solutions in the 
base-subsonic region are difficult to find, as discussed in ap- 
pendix^ In many cases, although shock solutions can be 
found, no solution of the correct periodicity exists. 

Many solutions show secondary peaks in the density. 
These occur especially near the ultraharmonic resonances, 
where higher harmonic terms in the solution become im- 
portant. If a secondary peak crosses the sonic line, a sec- 
ond shock will occur. It will therefore become impossible 
to locate a single-shock solution, since a singularity will be 
present at the second sonic point. It is possible to extend 
the method des cribed here to include a secondary shock 
JShu et alJl!973ft . but this has not been attempted in this 
paper. 



3.1.2 Results 

A 'standard model' is defined, with parameter values as 
given in tabled Density profiles for the standard model are 
shown in figure |S] for a range of radii. No profiles are shown 
for radii between 11 and 12.5 kpc, as no solutions exist. The 
figure shows the smooth transition toward the base-subsonic 
density profiles, in which the broad density peak at r\ = 
appears and the shock is small. This figure also indicates 
that the shock location moves to smaller values of r\ with 
increasing radius. At small radii, the shock occurs almost 
at the potential minimum, but as the radius increases and 
the shock becomes weaker, it moves upstream toward the 
potential maximum. 



3.1.1 Existence of a solution 

It is not always possible to find solutions containing a shock 
for any choice of parameters. The equations are based on the 
approximation that the gravitational field is in the r\ direc- 
tion (and hence that the shock is parallel to the equipo- 
tentials), which will hold for large or (equivalently) 

sin i <C 1, so that the spiral arms are tightly wound and $ (R) 
varies rapidly. In more open spirals the equations therefore 
break down. Also, the parameters are approximated as con- 
stant along a streamline, which will be a good approximation 



3.2 Numerical calculations 

The response of an isothermal gas disc can be computed 
with numerical hydrodynamical codes. In this section, re- 
sults obtained using two-dimensional Smoothed Particle Hy- 
drodynamics (2D SPH) and a two-dimensional Piecewise 
Parabolic Method (PPM) code are compared with the re- 
sults of the semi-analytical analysis. Numerical codes are not 
restricted by the approximations of the semi-analytical ap- 
proach, nor to results containing a single shock, and should 
therefore be able to capture more detailed structure. 
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Figure 5. Density profiles (7(77) for the standard model, at radii 
from R = 5 kpc to _R = 10.5 kpc (indicated) in steps of 0.5 kpc. 
The profile at R = 13 kpc is also shown. 

3.2.1 Two-dimensional Smoothed Particle Hydrodynamics 

Calculations with 2D SPH showed a good general agreement 
to the results of the semi-analytical analysis, but did not re- 
solve the sharp shocks accurately. A tendency for the SPH 
particles to 'clump' together unphysically was overcome by 
using a triquintic spline smoothing kernel, but t his results in 
decreased spatial resolution (see iGittina 12003 . for details). 
This problem may be alleviated by adjusting the time de- 
pendence of the smoothing l engths in the way described in 
lEnglmaier fc Gerhardl dl997tl . A typical example of results 
obtained with 2D SPH is shown in figure [(J 



3.2.2 Two-dimensional Piecewise Parabolic Method 

A 2D PPM code called cmhog2 was provided by P. Teuben. 
The code uses PPM in the Lagrangian remap for- 
malism, and includes isothermal hydrodynamics in po- 
lar coordinates. The cod e is b ased on that described 
in IPiner. Stone fc Teubenl lll995l) : for a description of 
the P iecewise Parabolic Method, see IColella fc Woodward! 

Results from calculations of the standard model and a 
model with m = 4 are shown in figure [7| Sharp shocks were 
formed as expected, and the strength of the shocks fades 
towards the corotation radius. Results from further calcu- 
lations, with variations of the strength F and the pattern 
speed S7p, were very similar. 

Results from the standard model are compared to semi- 
analytical results in figures El M and ESI at radii R = 10, 
7 and 13 kpc, respectively. Recall that a period of 27r in 
rj always corresponds to the passage from one arm to the 
next, regardless of the number of arms. The agreement is 
excellent in all three cases, and the narrow density peaks 
are reproduced by the grid code. The sharp shock front is 
spread over several zones, but it is reproduced much more 




R=10 kpc - 



-2 2 

V 

Figure 6. Density profile u(n) of the standard model at R = 10 
kpc as computed by the semi-analytical method (solid line), com- 
pared to results extracted from a 2D SPH calculation (circles). 




Figure 7. Surface density maps from calculations of the standard 
model (left) and with m = 4 (right) at t = 1 Gyr. The maximum 
radius is 20 kpc. 

accurately than in the SPH code. All the density profiles 
show some small oscillations in the density downstream of 
the shock. At R — 13 kpc, there is once again a discrepancy 
between the semi-analytical prediction and the calculation 
result, probably due to the proximity of the 6:1 ultrahar- 
monic resonance. 

The agreement with the semi-analytical profiles im- 
proves as the resolution of the grid-code is increased. Figure 
[SJshows the result obtained with 150 radial zones, compared 
to the full 600, for comparison. This trend held up to the 
limits of available resolution, and presumably, with more 
computing resources, ever more accurate agreement could 
be achieved. 

Figure 1111 shows a comparison at 10 kpc with the cal- 
culation using m = 4. The general shape of the profile is 
again well-represented, but in this case the position of the 
shock is slightly further upstream in the grid-code result 
than the semi-analytical prediction. This is a general result 
when m = 4. The reason for the difference is not clear. The 
number of grid zones covering one arm is half that in the 
m = 2 calculations, so doubling the resolution might remove 
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Figure 8. Density profile a(rj) of the standard model at R = 
10 kpc as computed by the semi-analytical method (solid line), 
compared to results from the grid-code calculation (circles). Also 
shown are the results from the same calculation with 1/4 the 
number of radial zones (crosses). 
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Figure 10. Comparison as in figurelSl at R = 13 kpc. This radius 
is in the base-subsonic region, and is close to the 6:1 resonance. 




Figure 9. Comparison as in figure El at R = 7 kpc. 

the offset. Such high resolution is impractical at time of writ- 
ing. However, when the angular resolution is halved instead, 
the predicted shock location is not obviously changed, which 
suggests that the effect is not resolution dependent. 

The shock location in the results corresponds very 
closely to the predicted location from the semi-analytical 
results in all cases where m = 2, at radii between the in- 
ner Lindblad resonance and the start of the base-subsonic 
region. Figure WI\ plots the measured location of the shock 
as a function of radius, from the calculation of the standard 




Figure 11. Comparison as in figure E| with m = 4. The result 
with half the angular resolution is also shown (crosses). 

model. The positions predicted from the semi-analytical the- 
ory are also shown, where they can be found. In figure Hoi 
the results for the model with m — 4 are given. When m — 2, 
there is excellent agreement between the two. At large radii, 
in the base-subsonic region, the shocks become so weak that 
the shock location in the grid-code results cannot be de- 
termined accurately. The comparison with m = 4 shows 
that the grid-code consistently places the shock further up- 
stream than the prediction of the semi-analytical theory, as 
discussed above. 

Similar comparisons were made from calculations using 
F = 0.03 and ttp = 19.5 km s _1 , and in both cases, very 
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Figure 12. Shock location vs. radius in the standard model, as 
measured from the grid-code calculation (points) and predicted 
from the semi-analytical results (crosses). 



Figure 14. Shock location with m = 2 measured from the grid- 
code calculation (points) over the range 5—15 kpc. Potential min- 
ima are marked as lines. Only the range < 9 < 7r is shown since 
the grid-code only calculates half the plane. 




close agreement was seen. Numerical results confirm the pre- 
dictions of the semi-analytical method in parameter regimes 
where such results are possible. Grid-code calculations were 
also carried out for parameters where no semi-analytical so- 
lutions exist, such as those with larger inclinations (i.e. more 
open arms) or lower sound speeds. In general, under these 
conditions, the gas failed to settle to a steady solution, and 
a much more complex structure was seen involving multiple 
shocks and significant secondary peaks in the density. These 
more complicated cases are left for further investigation. 



In figure ITU the shock location with m — 2 is plotted as 
lni? vs. 9 and compared to the potential minima, which are 
of course straight lines on such a plot. The shocks are near 
the minima at small radii, and move toward the maxima as 
corotation is approached, until they become too small for 
the grid-code to locate and the points scatter. The result 
with m = 4 is shown in figure H31 In this case, the shocks 
move from the maxima at small radii to cross the minima, 
and continue to move to the next maxima. The shocks thus 
complete a transition from one arm to the next. 
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4 SPIRAL ARM COMPONENTS AND THE 
OFFSET FUNCTION 

4.1 Shock locations 

The 'offset function' referred to in this paper, O(-R), is de- 
fined as 

= Tl(#shock — #min) 

where the potential minimum lies at angle 9 m i n and the 
shock is at angle # s hock- (The offset is defined to lie in the 
range — n < < 7r). For tightly wound spirals, in which 
the streamlines are almost circular, this is almost equivalent 
to r^hock — ??min. Since 7) is defined in such a way that the 
minima of the potential correspond to rj = (equation 0, 
rimxn = 0, and so 

= 7?shock- 

This very convenient result allows the offset to be read di- 
rectly from density profiles ct(t]), whether generated semi- 
analytically or in numerical calculations. 

The value of 0(-R) is difficult to predict without per- 
forming a full calculation. It will be expected to vary with 
the relative velocity of the gas and the potential, as it is the 
deceleration and acceleration of the gas that causes a shock 
to form in the first place. The offset should, therefore, vary 
systematically with radius in a galaxy, since the relative ve- 
locity R(£l — ttp) decreases with increasing radius towards 
corotation. 

The exist ence and varia tion of this offset have been 
noted before. iRoberts! |l969) noted, from calculations of 
semi-analytical solutions, that 'the shock lies just on the 
inner side of the background [i.e. potential] spira l arm', in- 
dicati ng that was generally small and negative. IShu et alJ 
il972T) report a value of = —72° for their calculat i ons o f 
the flow in the solar neig hbourhood. lYuan fc Grosb0ll Jl98lt) 
generated synthetic colour profiles across a spiral arm, based 
on semi-analytical solutions of the flow. They state that 'the 
shock . . . always occurs on the inner side of the potential 
minimum' (0 < 0), and for their calculations adopt a con- 
stant value = —30° (in their paper is given the symbol 
Ao). This result, that should be small and negative, is 
often tacitly assum ed to be correct i n disc ussions of spiral 
arms. For example, ISeiear fc Jamesl il998ll state that 'the 
fact that dust lanes [associated with shocks] appear on the 
trailing [i.e. upstream, inside corotation] edges of [potential] 
arms is itself evidence for the large-scale shock scenario'. 
The systematic variation of with radius, however, does 
not seem to have been discuss e d in th e literature. 

lElmegreen fc Thomassonl il993f) performed particle- 
based simulations of a galaxy, in which gas clouds and stars 
were represented by interacting particles. In the somewhat 
open (sini ~ 0.34) spiral arms that formed in the disc, they 
found that 'the location of the [gas] shock front changes 
from the inside of the [stellar] arm inside corotation to the 
outside of the arm outside corotation', that is, < inside 
corotation and > outside. This general trend is in agree- 
ment with the usual predictions discussed above. However, 
the shock front is very poorly resolved in their results, and 
a d etailed examination of would n ot be possible. 

iKikuch i. Korchagin fc Mivamal l)l997|) discuss angular 
separations between different arm components arising in an 
unstable galactic disc. Their analysis differs from ours in 



that they constrain all components to share a common ve- 
locity field, thus negating the possibility of shocks in the gas 
phase. They nevertheless find that the ordering of the gas 
and stellar arms reverses at corotation, and suggest the util- 
ity of this result in constraining the location of corotation 
in observed spirals. 

4.2 Observational characteristics of spiral arms 

In discussing the observational characteristics of a model 
spiral galaxy, it is very important to be precise about what 
is meant by an 'arm'. Since a spiral arm will affect different 
components of the disc in different ways, arms traced by 
one phenomenon may appear very different to arms traced 
by another. 

In general, one expects to see three spiral arm patterns. 
The first is that traced by the potential minima, or the max- 
imum surface density. In fact, these two definitions will not 
coincide exactly, since the potential is always more curved in 
R than the underlying density. The effect is that the poten- 
tial minima are slightly more tightly wound than the peaks 
of the surface density, but we will disregard this difference 
here. These arms, created by the mass variations in the old 
stellar population, are expected to have a smoothly varying 
structure. The second set of arms is traced by the shock 
front of the gas (where it exists). In strong shocks, this will 
also represent the maximum gas density. Since dust is gen- 
erally taken to be a tracer of shocked gas, these arms will 
also be traced by dust lanes. The third set of arms is the re- 
gion of enhanced star formation. The existence of this third 
arm, linked to the other two, relies on the assumption that 
the shocked gas promotes star formation in the ISM. There 
will then be some finite 'onset' time, during which the early 
stages of star formation take place. This is followed by an in- 
creased density of young stars, Hn regions and so on, which 
will trace the region of star formation. This third set of arms 
would, therefore, be expected to occur downstream of the 
second. These three arms are labelled, in this paper, the P- 
arm (for Potential), the D-arm (for Dust) and the SF-arm 
(for Star Formation). 

With this definition, refers to the offset between the 
P-arm and the D-arm, irrespective of the existence or loca- 
tion of any SF-arm. If 0(i?) is a constant, then the D-arms 
will be identical to the P-arms, but rotated by an angle 
0/m, and the two will therefore have identical pitch angles. 
If, however, becomes more negative with increasing ra- 
dius (e.g. the shock moves further upstream of the potential 
minimum at larger radii), then the D-arm should be more 
tightly wound than the P-arm. 

The relative pitch angles of the D-arm and SF-arm are 
also likely to be different. If the assumption of a constant 
time offset from the D-arm to the SF-arm (associated with 
the onset of star formation) is made, then the angular offset 
will be proportional to the relative angular velocity of the 
gas SI — Qp. Since this decreases to zero as radii approach 
corotation, the offset will also decrease to zero. The SF-arm 
would therefore be expected to be more tightly wound than 
the D-arm. 

Figure [Tol is a diagram of the three arms, in the partic- 
ular case that becomes more negative with radius and the 
star formation onset time is constant. Under both of these 
circumstances, the pitch angles of the arms are expected 
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SF-arm (star formation) 




Figure 16. Diagram of the three types of spiral arm and their 
separations. In this diagram © becomes more negative with ra- 
dius, so that the D-arm is more tightly wound than the P-arm. 
The star formation time delay is constant, but since the relative 
velocity of the gas to the P-arm is smaller at larger radii, the an- 
gular delay is smaller. Consequently the SF-arm is more tightly 
wound than the D-arm. 



to follow ip > in > isF- IVisserl <ll980all makes the point 
that a constant O is required for the P-arms and D-arms to 
have the same pitch angle, and comments that in M81 the 
tendency is for both the 'gas and dust' (D-arms) and 'H n 
regions and young stars' (SF-arms) to be tighter than the 
potential arms. 

The three types of spiral arm are observed in different 
ways. Observing the P-arm requires observations at longer 
wavelengths, that isolate the older stellar population in the 
disc. These arms would probably, therefore, be observed as 
peaks in the K or I bands. The D-arms, defined by the lo- 
cations of shocks, should be traced by dust lanes or radio 
emission ridges. If these arms correspond to the peaks of 
the gas density, then they should also be traced by peaks of 
CO emission (and other tracers of dense molecular gas). Fi- 
nally, the SF-arms can be traced by bluer light, such as the 
B band, or by Ha emission from Hn regions. Determining 
the location and extent of the SF-arms is hampered by the 
effects of dust reddening. 

Generally, the strongest CO emission is seen to coincide 
with t he dust lanes an d nonthermal ridges, in M51 for ex- 
ample (iLo et al.lll987F) . an d all of these are displaced from 
the p eak of Ha emission jTilanus et al.lll988l : iNakai et alJ 
1994). The offset from CO peaks to Hn regions is seen 
in m ost grand design spirals, such as M31 ijlchikawa et alJ 
Il985l) . More recent observations have detected a sequence of 
12 CO to 13 CO to Ha across an arm in M51 llTosaki et alJ 
2002), supporting the idea of a gradual collapse of molecular 
clouds to higher densities in the onset of star formation . One 
exception is the observation of lLord fc Kennevl (Il99 lTl that 
peaks of CO emission are offset from the dust lanes in an 
arm of M83, which they attribute to dense clouds penetrat- 
ing through the shock layer while a more diffuse component 
is compressed to a shock. 

In principle, therefore, the positions of the D-arm and 
P-arm can be observed. This would then allow a direct mea- 
surement of Q(R). The resulting data will enable constraints 



to be placed on the parameters of the galaxy's spiral struc- 
ture, as discussed in the next section. 



5 LOCATING THE COROTATION RADIUS IN 
GALAXIES 

The determination of the corotation radius (or, equivalently, 
the pattern speed) in spiral galaxies has been attempted by 
various methods for various galaxies. The assumption that 
a pattern speed exists at all implies that a rigidly rotating 
spiral potential is present. 

The simplest way to estimate the corotation radius is 
to associate it with the maximum extent of observed spiral 
arms. The idea that star formation is promoted by the shock 
compression means that at corotation (where there can be 
no shock compression), star formation should be absent. The 
radius of the outermost H II region has there fore been u sed 
as an estimate of the corotation radius (e.g. iRotd Il975l in 
M81). This method is not reliable, since longer exposures 
tend to reveal further features at larger radii in galaxies, 
and H II regions may exist outside corotation. An alterna- 
tive, but similar, approach is to associate the inner edge of 
observed spiral features with the inner Lindblad resonance, 
on the basis that spiral waves are not expected to extend 
beyond the Lindblad resonances. However, several investiga- 
tions have discussed the possibility of spiral patterns extend- 
ing within the inner Lindblad resonance jElmeg reen et alJ 
Il99ci lEnglmaier fc Shlosmanll200oT: iMartini et al.ll2003ft . "so 
this method is unlikely to be accurate without detailed mod- 
elling. More sophisticated analyses associate various obser- 
vational features with a range of resonances in a galaxy 
to constrain the corotation radiu s jElmegreen fc, Elmegreenl 
ll99CtlElmegreen. Elmegreen fc Montoiegrol^9^ 

Alternatively, if the velocity curve (and preferably other 
parameters such as velocity dispersions and mass models) of 
a galaxy are known, then a corresponding spiral structure 
can be simulated, based on a specified pattern speed and 
the spiral density wave theory. The pattern speed can then 
be adjust ed unti l the bes t fit is achieved. This approach was 
taken by IVisserl Jl980bh in an analysis of M81, and more 
recently for several galaxies jKranz. Slvz fc Rbjfeooj) . 

Another approach makes use of the predicted radial 
variations in the arm amplitude. These are associated with 
interference between outward and inward travelling spi- 
ral w aves, and fitting the oscillations can locate corota- 
tion (Elmcgrccn, Scidcn & Elmcgrccn 1989). The 'modal' 
approach (Be rtin et alJ 11989) involves examining the full 
range of normal modes of the galactic disc and their rel- 
ative amplitudes. The technique has been applied to M81 
jLowe et alJll994h . Other methods include: matching ra- 
dial oscillati ons in the velocity curve to predicted ve- 
locity fields dYuanl Il969l) : the 'geometric phase' method 
jCanzian fc Alien! Il99/tl . which associates corotation with 
a shift from singly to triply symmetric structure in the ve- 
locity field; and the integrated continuity equation method 
( Wcstpfahl 1998), in which asymmetries in the observed disc 
are exploited to constrain Qp. 

Of course, values determined by different methods 
rarely agree exactly; table|5|lists pattern speeds derived for 
the spiral pattern in M81 from a variety of sources in the 



10 D. M. Gittins and C. J. Clarke 

Table 2. Some pattern speeds of M81 reported in the literature. 



Qp I km s 1 kpc 



Source 



17.4 

26.4 

20 

18.3 

18 

20 

17 ±2 

26 

24 

23.4 ±2.3 



Shu. Stachnik fc Yost (19711 
Roberts. Roberts fc Shu (19751 
Rots (19751 

Gottesman & Wcliachcw (1975) 
Visser (1980a 1 
^errrian^^isjrurov 
Sakhibov fc Smirnov ( 19871 



Elmcg reen et al. (19891 
Lowe et al. (19941 



Wcstpfahl (1998) 



literature. This indicates a spread of some 25% in one of the 
most studied spiral galaxies. 

A commonly suggested method of locating corotation 
arises from examination of the separate arm components. If 
star formation occurs at some time offset after the passage 
of a spiral shock, then the relative location of the D-arm 
and the SF-arm should be reversed inside and outside coro- 
tation. The radius at which these arms cross should, there- 
fore, correspond to coro tation. This point is mentioned by 
iGrosbal fc Patsis! 1^998), but they failed to observe such a 
cros sing. 

lDixonl ( ll97lD proposed that observations in the B and V 
bands could allow the variation of stellar age across an arm 
to be determined, and that corotation should correspond to 
the radius at which this trend reverses angular direction. 
IPuerari fc Dottoril il997tl attempted to locate corotation in 
two galaxies by this method, performing a Fourier analy- 
sis on observations in B and I bands. Observational studies 
of the separations of arm components usually focus on the 
SF-arm a nd the P-arm, observing in (for example) B and 



K bands. iBeckman fc Cepal lll990D observed in the B and 
I bands, finding an offs et in one galaxy, but too flocculent 
a structure in another. iGrosbol fc Patsis! jl998f) measured 
radial variations in the location of SF-arms and P-arms, 
but conclude that uncertainties in the absolute phase off- 
set between them pre vent corotation from being located. 
ISeiear fc Jamesl (Il998l) also examined arms in the K and 
B bands, observing that the SF-arms appear on the con- 
vex side of the P-arms as expected, but that the P-arms are 
more tightly wound than the SF-arms, which is the opposite 
of the expected result. 

All of these investigations are based on observations of 
the SF-arms. This will always be inherently difficult, since 
the existence, location and shape of the SF-arm are neither 
expected, nor observed, to be very regular. As discussed 
in the introduction, galaxies with regular structure in P- 
arms will often be classified as flocculent in optical obser- 
vations, meaning t hat they do not p ossess regular SF-arms. 
As pointed out bv lElmegreenl il979Tl . one would expect the 
stimulation of star formation in shock fronts to be followed 
by a 'shuffling' of material, including feedback from star for- 
mation, random structure from the condensation of dense 
clouds and so forth. The appearance of any SF-arm would 
therefore be much more irregular, and in the right circum- 
stances will appear entirely flocculent. The observation of 
these arms is further hampered by the effects of dust obscu- 
ration, which can make it difficult to locate their full extent. 
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Figure 17. Plots of &(R) with various parameters. Lines marked 
with crosses correspond to (bottom to top) F = 0.03, standard 
model, i = 0.15, a = 5 km s _1 , m = 4. Line with circles corre- 
sponds to ftp = 19.5 km s — 1 kpc -1 . 



Furthermore, the unknown quantity of the 'onset time' for 
star formation may in reality be a range of times, leading to 
star formation occurring at a range of offsets from the shock. 
The prediction of the large scale shock-induced star forma- 
tion scenario would be, therefore, that regular P-arms should 
produce regular D-arms, and any subsequent SF-arms will 
be irregular and poorly localised. 

Observations of the P-arms and D-arms, and the offset 
O between them, avoid all such issues and would be ex- 
pected, in galaxies with a regular spiral potential, to result 
in a much clearer result. In the next section, the way in which 
corotation might be constrained from such measurements is 
described. 



5.1 Constraining corotation from the 
shock—potential offset 

As has already been seen in figure 1121 there is a general 
trend for O to vary in two-armed spirals in a specific way. At 
small radii, O approaches zero, but as the radius approaches 
corotation, O moves toward — n. Even with m — 4, although 
O does not go to zero at small radii, it still moves toward — n 
with increasing R. The extrapolation of the curve of O, in 
all these cases, roughly indicates the location of corotation. 

The form of Q(R) does, of course, depend on the param- 
eters of the spiral potential. However, the extrapolation of O 
to find corotation stays approximately unchanged under all 
such variations. Figure [T71 shows the curve for the standard 
model compared to the alterations i — 0.15, F = 0.03, a = 5 
km s , m — 4 and Qp = 19.5 km s _1 kpc" 1 . Corotation is 
located at 17 kpc in all of these models except the last, in 
which it is at 11.3 kpc. The shapes of all the curves indicate 
corotation to be in the range 15-20 kpc, except the curve 
at higher fip which correctly indicates a smaller radius of 
corotation. 
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Figure 18. Velocity curve used in lShu et alj ll97i 



Variations in the velocity curve will also affect the shape 
of 0(7?). In figure Hoi the resulting plots are compared using 
a flat velocity curve, a rising velocity curve in which e'j -1 = 10 
kpc (rather than 1.5 kpc), and the velocity curve used in 
IShu et all <ll973f) . which is peaked at around 8 kpc, and 
shown in figure [TH1 The values corresponding to fip = 19.5 
km s kpc -1 are also shown. Since each of these veloc- 
ity curves places corotation at a different radius, the values 
are plotted against -R/iicorotation in this figure. Again, in all 
cases the approximate position of corotation is indicated by 
the extrapolation of the curve. 

Probably the most difficult parameter to constrain in a 
spiral galaxy is the relative strength F of the spiral pertur- 
bations. In the standard model, this takes the value 5% at 
R = 8.5 kpc, and varies as Rexp(—e s R). In figure |2T)1 the 
shock locations obtained using a constant strength F(R), 
at values 3%, 5% and 10%, are shown. This figure indi- 
cates that the spiral strength introduces the greatest un- 
certainty into the location of corotation. Generally, weaker 
spirals will form shocks further upstream at smaller radii 
than stronger spirals. If the strength of the spiral perturba- 
tion varied systematically across a galactic disc, for example 
increasing from 3% to 10% (or vice versa), the extrapolation 
of © could indicate an inaccurate value for corotation. 
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Figure 19. Plots of S(R), plotted against radius relative to coro- 
tation, for various parameters. Lines marked with crosses corre- 
spond to (left to right) e^ 1 = 10 kpc, standard model, velocity 
curve from figure 1181 Line marked with circles corresponds to 
ftp = 19.5 km s" 1 kpc" 1 . 
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5.2 Example results and accuracy 

To assess the approximate accuracy of the method, we con- 
sider in this section how the semi-analytical results in figures 
I17lll9l and l20l mav be used, in conjunction with observational 
data, to constrain the location of corotation in spiral galax- 
ies. We therefore begin with 10 sets of data Q(R) for m — 2. 
The general idea is to extrapolate the downward-curving 
trend to the point © = — n, which forms the estimate for 
corotation. 

The radial range over which we have managed to obtain 
semi-analytic solutions depends on the model parameters 
(e.g. arm strength, rotation curve), so that these parameters 



Figure 20. Plots of @(R) for models using a constant F of (bot- 
tom to top) 0.03, 0.05 and 0.1. 

affect the range in Q(R) for which we have solutions in each 
case. Estimates will be more accurate where a larger A© 
is available. We wish to assess the approximate accuracy of 
this method as a function of AO. 

To this end, we analysed the results at values of AO in 
steps of 7r/8 up to the largest value covered by the data. For 
each value of AO, only those curves were included whose 
results cover at least that range, and only the points (at 
small radii) within that range were used. Thus, at small 
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A0/(tt/8) 



Figure 21. The extrapolation of the results with a = 5 km s _1 
at two values of A©. At A© = 7r/8, only the first five points 
constitute the data (solid circles). The resulting quadratic fit has 
the wrong curvature, so a linear fit is used (dashed line). At A© = 
it/2, all the data are included, and the quadratic fit estimates the 
corotation radius at 18.6 kpc. 



AO, all the curves can be extrapolated, but only the first 
part of each will be used. At large AO, only a few of the 
curves have sufficient data to produce a result, but they will 
be more accurate. 

Any set of data clearly needs to show a decrease in 
with R as a minimum requirement for this method to 
be applied. The simplest extrapolation is to apply a least- 
squares quadratic fit to the points. In some cases, the data 
curve toward positive O over some sections, in which case 
a quadratic fit does not extrapolate to O = — tt. To avoid 
this problem, if the coefficient of R? in the quadratic fit 
is positive, we use a linear fit instead. Once a fit has been 
chosen, the corotation radius is estimated as the intersection 
of the fit with O = — n. Figure 1211 illustrates some fits for 
the data with a — 5 km s _ . 

The scatter in the resulting predicted values for coro- 
tation over the 10 sets of data gives an indication of the 
accuracy of this method. The results are plotted in figure 
1221 for all values of AO. At an optimal value of AO ~ 3n/8 
to 47r/8, most of the data sets cover a sufficient range to be 
included, and the spread is some 25%. We estimate, there- 
fore, that this method should be accurate to around 25% if 
a sufficiently large AO is available in the observed data, and 
that AO >7r/4 is a reasonable guide to this requirement. 



6 CONCLUSIONS 

In this paper, we re-examine the response of an isother- 
mal gas disc to a model spiral potential, leading to a large- 
scale flow containing a single shock. The established semi- 
analytical method was compared to results from numerical 
techniques, with good agreement within the limitations of 
each method. We found that: 



Figure 22. The ratio of predicted to actual corotation radius for 
all the data sets, as a function of A©. Solid lines with crosses 
correspond to (bottom to top at left end) F = 0.03, standard 
model, a = 5 km s — 1 , i = 0.15. The solid line with circles cor- 
responds to Qp = 19.5 km s — 1 kpc -1 . Dashed lines are models 
with constant F with values (bottom to top at right end) 3%, 5%, 
10%. Dotted lines are alternative velocity curves (lower: t~ l = 10 
kpc, upper: figurc llSl . Solid circles indicate results where a linear 
(rather than quadratic) extrapolation was used. 



(i) Results from two-dimensional SPH fail to resolve the 
shock sufficiently sharply to determine its location. More 
work on the details of the SPH implementation may improve 
results. 

(ii) Much better results were obtained with the PPMLR 
code CMHOG2. 

(iii) In parameter regimes where a semi-analytical solu- 
tion cannot be found, in general the numerical results failed 
to show a simple single-shock solution. 

We discuss the appearance of spiral arms in galaxies, 
and the existence of three separate arms, the P-arm (in the 
stellar disc, approximately tracing the gravitational Poten- 
tial), the D-arm (the shock as traced by Dust) and the SF- 
arm (an arm traced by enhanced Star Formation). Even in 
galaxies classified as optically flocculent, a regular structure 
in the P-arms and D-arms should be common. The angular 
offset 0(i?) between the D-arm and P-arm is proposed as 
an observable function, avoiding all the difficulties inherent 
in tracing the shape and location of the SF-arm. 

Finally, we show how this function could be extrapo- 
lated to provide an estimate of the corotation radius, and 
assess the accuracy of the method as a function of the range 
AO available in the observations. We estimate that given 
AO>7r/4, the method is accurate to around 25%. 
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APPENDIX A: METHOD FOR LOCATING 
SEMI-ANALYTICAL SOLUTIONS 

In this appendix, the method for locating solutions in the 
semi-analytical approach is described for conven ience. The 
metho d is almost the sam e as that described in IShu et alJ 
Jl972l) and lRobertsl Jl969h . 

The differential equations describing the response of the 
gas, along a streamline with average radius Ro, are as fol- 
lows: 

du„i n 2itf i — FRqQ sin r\ 
— — = U(u v o +u v i)- ■ ^ j- 

ar] (u v0 + u vl ) 2 - a 2 

gMgl _ y 

drj u n0 + u v i 

where 

u v o = Ro(Q. — f2p) sini 
u^o — Ro(Sl — f2p) cosi 

u = ^R n 

m 

_ sini d /p3n\ _ sini R k 2 
~ m diT ' ~ m 2Q 

The resulting gas density a = ctq + o\ relative to the 
unperturbed density ao(R) is given by 

(<T + (Ti)(m^o + U n i) = a u v0 . 

It is convenient to transform these equations into di- 
mensionless form with (u v i, u^i) replaced by the dimension- 
less (w, v ) defined by 

U — , 
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by first defining the dimensionless parameters 
FRofl 



f 



2V 

—u v o 



2UV 

so that the final dimensionless equations to be solved are as 
follows: 



du , . v ■ 
— = (u-v)-— 
Or) (it 



dr) 



/sin?? 



(Al) 



(A2) 



For low spiral strength F or sufficently high or low 
sound speed a, solutions will be entirely supersonic or en- 
tirely subsonic, and will not contain a shock. For parameters 
where u v crosses the value a, a shock must exist. The method 
for locating such a solution containing a single shock (if it 
exists) is as follows. 

As can be immediately seen from equations IAll and lA2l 
singularities occur at three points: u = v, u — v — y/x and 
u — v-\-\fx. These correspond, respectively, to u v = 0, u v = 
—a and u v = a. For a solution flowing through the shock, 
the first two conditions should never occur, since gas flow 
through a trailing spiral (inside corotation) should always be 
moving outward in 77. The third singularity corresponds to 
the sonic point, at which the gas undergoes the transition to 
supersonic flow after the shock (the gas returns to subsonic 
flow in the shock itself). The value of r\ at which this sonic 
point occurs is denoted r/sp. The singularity at this point 
is removed as long as v(r)sp) = /sin^sp. This requirement 
fixes the boundary conditions; given any guessed value for 
77SP, the starting values of u and v at this point must be 
u = v + \fx and v = / sin 77SP . 

Solutions are therefore calculated starting at r)sp, and 
proceeding forward in 77 (the subsonic branch) and back- 
ward in r\ (the supersonic branch). The determination of 
the shock location and the procedure for adjusting jjsp to 
find a solution are described below. 

The presence of the singulari ty at r)sp can cre ate dif- 
ficulties for numerical integration. IShu et al] il973l) used a 
series expansion in 77 around ?isp to overcome this problem. 
An alternative is to rewrite the equations in new variables, 
which hide the singularity, as follows (J. Gair, private com- 
munication). Start by defining 

u — v + \fx + u 

v — f sin 7) + v 

and then define 



The differential equations now become 

dw _ 2 (Vi ± Vwj v 

9v 2^/x±Vw 

dv v 

— — = — 1 — f cos ri j= 

orj ■ s [x± Vw 

where the positive root of w is taken where u > (i.e. along 



the supersonic branch) and the negative root is taken where 
u < (the subsonic branch). In this form, the equations 
can be integrated along each branch separately, starting at 
77SP, and w and v transformed back to find u and v. The 
singularity is hidden and replaced by a forbidden region of 
phase space, namely the requirement that w > 0. This will 
be satisfied so long as the correct boundary condition is used 
(v(r]sp) = w(??sp) = 0). 

The procedure for locating the shock is as follows. The 
flow is integrated from /?sp along the supersonic and sub- 
sonic branches, to give u aU p(r]) and u S ub(r])- Generally, the 
supersonic branch will eventually return to the singularity 
at u = a and cannot be further integrated, while the sub- 
sonic branch either tends asymptotically to u — 0, or also 
reaches u = a. If a consistent solution exists, then the shock 
will terminate both branches before they reach singularities. 

For every point along the supersonic branch u S u P , the 
isothermal jump conditions can be used to calculate the 
post-shock velocity Mahock that the flow would have if the 
shock were to occur at that particular value of 77. Since 
the shock is perpendicular to the rj direction and the gas 
is isothermal, the jump condition is simply 

— r, 2 

T^sup^shock — & 

and so the post-shock branch is defined (in terms of the 
dimensionless parameters) by 

Ushock^) = V + 



Through the shock, v should be continuous, but u should 
drop from M sup to M a hock- In a consistent solution, the shock 
joins the supersonic branch to the subsonic branch at the 
same value of v. Therefore, the only consistent location for 
the shock exists where the post-shock branch crosses the 
subsonic branch in the (v,u) plane. This process is shown 
visually in figure IA1I The starting value of ?7sp must be 
adjusted until a shock location appears, at a point on each 
branch before they reach a singularity. For some parameters, 
no value of ?7sp will produce a shock. 

The location of the crossing point is found numerically 
by minimising u sub (rj sub ) - u shock (r] SU p) with respect to 77 sub 
and 77su P , which are the values of 77 on the subsonic and 
supersonic branch (respectively) at which the shock occurs. 
The period in 77 of the solution will therefore be rj SU p—rj sub . In 
general, for a given value of the sonic point ?7sp, this period 
will not be 2n. The starting value of ?7sp must therefore 
be adjusted, repeating the shock location procedure at each 
value, until a solution is found with a period of 2ir. If such a 
value can be found, it represents the only physical solution 
to the flow equations, containing one isothermal shock. The 
density a(rj) can be calculated from the values of u. 

Solutions in the base-subsonic region are particularly 
hard to find. This is because the integrated branches depend 
very sensitively on the value of ?7sp. In some cases, the range 
of values for which a solution containing a shock exists can 
have a width of less than 7r/1000, and so can be difficult to 
locate. Outside this range, the calculated branches run into 
singularities too rapidly. 
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Figure Al. The process of locating the shock. Integration starts 
at the sonic point (labelled »7sp) along the supersonic (dotted) and 
subsonic (short dashed) branches. The post-shock branch (solid) 
is calculated for each point along the supersonic branch using the 
jump conditions. It crosses the subsonic branch at the appropri- 
ate location for the shock (circle). The shock, indicated by the 
long dashed line, closes the velocity curve. Flow proceeds in an 
anticlockwise direction. The value of r/sp should now be adjusted 
until the period in r\ of the solution reaches 2n. The unperturbed 
flow corresponds to u = v = (cross). 



